#### prepare functions ####
# function for equivalene testing
equivarence.test <- function(first.group,  # data frame of the first group
                             second.group,  # data frame of the second group
                             covariate,  # covariates' names
                             epsilon  # equivalence limit
) {
  x.t <- first.group[, covariate]  # covariate in the first group
  x.c <- second.group[, covariate]  # covariate in the second group
  n.t <- sum(! is.na(x.t))  # omitting missing values
  n.c <- sum(! is.na(x.c))
  mean.t <- mean(x.t, na.rm = TRUE)  # mean in the treatment group
  mean.c <- mean(x.c, na.rm = TRUE)  # variance in the control group
  var.t <- var(x.t, na.rm = TRUE)  # mean in the treatment group
  var.c <- var(x.c, na.rm = TRUE)  # variance in the control group
  pooled.sd <- sqrt(((n.t - 1) * var.t + (n.c - 1) * var.c) / (n.t + n.c - 2))  # pooled standard deviation
  t.stat <- sqrt(n.t * n.c * (n.t + n.c - 2) / (n.t + n.c)) * (mean.t - mean.c) / 
    sqrt(sum((x.t - mean.t) ^ 2, na.rm = TRUE) + sum((x.c - mean.c) ^ 2, na.rm = TRUE))  # test statistic
  ncp <- (n.t * n.c * epsilon ^ 2) / (n.t + n.c)  # noncentrality parameter
  std.diff <- (mean.t - mean.c) / pooled.sd  # standardized difference in means
  p.values <- pf(abs(t.stat) ^ 2, 1, n.t + n.c - 2, ncp)  # p-value in equivalence testing
  c(std.diff, p.values)
}

# wrapper function for equivalence tests
EQT.wrapper.conjoint <- function(data.tested,  # data frame containing variables
                                 treatment,  # treatment variable's name
                                 covariates,  # covariates' names
                                 epsilon  # equivalence limit
) {
  treatment.vec <- data.tested[, treatment]  # vector of the treatment variable
  agree.group <- data.tested[treatment.vec == "Agree", ]  # agree condition's data frame
  neutral.group <- data.tested[treatment.vec == "Neither", ]  # neither condition's data frame
  disagree.group <- data.tested[treatment.vec == "Disagree", ]  # disagree condition's data frame
  results <- matrix(NA, length(covariates), 6)
  for (i in 1:length(covariates)) {
    # agree v. neutral
    results[i, 1:2] <- equivarence.test(first.group = agree.group, second.group = neutral.group, 
                                        covariate = covariates[i], epsilon = epsilon)
    # agree v. disagree
    results[i, 3:4] <- equivarence.test(first.group = agree.group, second.group = disagree.group, 
                                        covariate = covariates[i], epsilon = epsilon)
    # neutral v. disagree
    results[i, 5:6] <- equivarence.test(first.group = neutral.group, second.group = disagree.group, 
                                        covariate = covariates[i], epsilon = epsilon)
  }
  # p-value adjustment by the Benjamini-Hochberg method
  for (i in c(2, 4, 6)) {
    results[, i] <- p.adjust(results[, i], method = "BH")
  }
  rownames(results) <- covariates
  colnames(results) <- c("std.dif.1", "p.value.1", "std.dif.2", 
                         "p.value.2", "std.dif.3", "p.value.3")
  results
}

#### load data ####
task.data <- read.csv("task_data.csv")

# reorder the levels of attribute variables
position.label <- c("Agree", "Neither", "Disagree")
task.data$a.article9 <- factor(task.data$a.article9, levels = position.label)
task.data$b.defense <- factor(task.data$b.defense, levels = position.label)
task.data$c.revisionism <- factor(task.data$c.revisionism, levels = position.label)
task.data$d.women <- factor(task.data$d.women, levels = position.label)
task.data$e.gay <- factor(task.data$e.gay, levels = position.label)
task.data$f.immigrant <- factor(task.data$f.immigrant, levels = position.label)
task.data$g.growth <- factor(task.data$g.growth, levels = position.label)
task.data$h.tax <- factor(task.data$h.tax, levels = position.label)

#### analysis ####
# Table A.3(a)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "a.article9", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(b)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "b.defense", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(c)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "c.revisionism", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(d)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "d.women", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(e)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "e.gay", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(f)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "f.immigrant", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(g)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "g.growth", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)

# Table A.3(h)
round(EQT.wrapper.conjoint(data.tested = task.data, treatment = "h.tax", 
                           covariates = c("gender", "age", "low.edu", 
                                          "middle.edu", "high.edu", 
                                          "DID.ratio", "ideology", 
                                          "article9.attitude", "defense.attitude", 
                                          "revisionism.attitude", "women.attitude", 
                                          "gay.attitude", "immigrant.attitude", 
                                          "growth.attitude", "tax.attitude", 
                                          "knowledge"), 
                           epsilon = 0.36), 3)